Environmental  Laboratory  erdc/el tr-09-4 


US  Army  Corps 
of  Engineers® 

Engineer  Research  and 
Development  Center 


System-Wide  Water  Resources  Program 

Development  of  a  River  and  Stream  Water 
Quality  Module 

Billy  E.  Johnson  and  Zhonglong  Zhang  June  2009 


Approved  for  public  release;  distribution  is  unlimited. 


System-Wide  Water  Resources  Program 


ERDC/EL  TR-09-4 
June  2009 


Development  of  a  River  and  Stream  Water 
Quality  Module 


Billy  E.  Johnson 

Environmental  Laboratory 

U.S.  Army  Engineer  Research  and  Development  Center 
3909  Hails  Ferry  Road 
Vicksburg,  MS  39180-6199 

Zhonglong  Zhang 

SecPro  Inc. 

3909  Halls  Ferry  Road 
Vicksburg,  MS  39180-6199 


Final  report 

Approved  for  public  release;  distribution  is  unlimited. 


Prepared  for  Fleadquarters,  U.S.  Army  Corps  of  Engineers 
Washington,  DC  20314-1000 


ER  DC/EL  TR-09-4 


ii 


Abstract:  This  report  describes  the  in-stream  water  quality  module 
within  the  Nutrient  Sub-Model  (NSM),  Version  1.2.  The  in-stream  water 
quality  module  includes  the  major  processes  for  nitrogen,  phosphorus, 
carbon,  and  oxygen  cycling  in  a  stream.  Kinetic  process  equations  are 
presented  for  each  state  variable  for  water  quality  modeling  in  streams. 
These  equations  were  taken  largely  from  models  such  as  QUAL2K,  ICM, 
etc.  and  summarize  past  development  efforts.  As  research  continues,  it  is 
anticipated  that  improved  process  descriptions  will  be  developed  and  will 
be  integrated  into  the  NSM. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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Preface 

This  report  summarizes  the  efforts  undertaken  in  developing  a  nutrient 
sub-model  for  linkage  with  a  variety  of  hydraulic  and  hydrologic  modeling 
systems  since  the  release  of  SWWRP-NSM  version  1.1.  This  development 
effort  was  performed  by  the  U.S.  Army  Engineer  Research  and  Develop¬ 
ment  Center  (ERDC),  Vicksburg,  Mississippi.  Funding  was  provided  under 
the  System-Wide  Water  Resources  Program  (SWWRP).  Appreciation  is 
extended  to  all  those  who  assisted  in  the  formulations  and  review  of  the 
process  descriptions  implemented  within  the  Nutrient  Sub-Model  (NSM). 

Principal  Investigators  for  this  study  were  Dr.  Billy  E.  Johnson  of  the 
Water  Quality  and  Contaminant  Modeling  Branch  (WQCMB),  Environ¬ 
mental  Laboratory  (EL)  and  Dr.  Zhonglong  Zhang  of  SpecPro  Inc. 

Dr.  Zhang  was  funded  as  an  onsite  contractor  under  Task  Order  Contract 
W912HZ-05-D-0011  on  Civil  Delivery  Order  No.  15.  Dr.  Johnson  con¬ 
ducted  his  portion  of  the  study  under  the  general  supervision  of  Dr.  Barry 
W.  Bunch,  Chief  of  WQCMB,  and  under  the  general  supervision  of 
Dr.  Richard  E.  Price,  Chief,  Environmental  Processes  and  Engineering 
Division,  and  Dr.  Beth  C.  Fleming,  Director  of  EL. 

COL.  Gary  E.  Johnston  was  Commander  and  Executive  Director  of  ERDC. 
Dr.  James  R.  Houston  was  Director. 
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1  Introduction 

Within  the  System-Wide  Water  Resources  Research  Program  (SWWRP), 
multiple  riverine,  estuarine,  watershed,  and  subsurface  flow  models  are 
being  modified  to  address  issues  of  environmental  concern.  Several 
integration  approaches  are  either  ongoing  or  have  been  proposed  to 
accomplish  this  task.  To  have  a  full  system-wide  water  quality  and 
contaminant  capability  in  SWWRP,  the  various  hydrologic  and  hydro- 
dynamic  engines  must  utilize  a  common  water  quality  and  contaminant 
approach  to  prevent  the  arbitrary  partitioning  of  constituents.  The  goal  of 
this  development  effort  is  to  upgrade  existing  hydrologic  and  hydraulic 
models  (i.e.,  water  engines)  using  a  common  water  quality  approach  in 
order  to  facilitate  their  linkage  and  application  on  a  system-wide  basis. 

In  keeping  with  a  common  water  quality  approach  to  model  development, 
a  library  of  water  quality  kinetic  modules  are  being  developed  such  that 
they  can  be  integrated  with  a  variety  of  water  transport  engines.  Modules 
in  the  library  of  algorithms  have  the  following  characteristics: 

•  multi-species,  multi-phase,  and  multi-reaction  system 

•  fast  (equilibrium-based)  and  slow  (non-equilibrium-based  or  rate- 
based)  reactions 

•  easily  extensible  to  new  reaction  pathways 

•  include  both  common  nutrient  and  contaminant  packages  as  well  as 
geochemistry 

•  simple,  well-defined  data  interface  and  calling  procedure 

The  water  quality  modules  are  developed  such  that  they  are  data  structure 
independent  thus  facilitating  their  integration  into  a  wide  range  of 
modeling  systems. 

The  Nutrient  Sub-Model  (NSM)  has  been  developed  as  a  library  of  water 
quality  kinetic  modules.  NSM  considers  a  detailed  N  and  P  cycling  and 
computes  nutrient  kinetic  fluxes  for  nitrogen,  phosphorus,  and  other 
species  in  watersheds  to  water  bodies  (Johnson  and  Gerald  2006;  Johnson 
et  al.  2007).  Modeling  of  nutrients  in  previously  released  versions  of  NSM 
consists  of  three  distinct  parts:  1)  Simulation  of  the  major  N  and  P  cycling 
processes  in  the  soil;  2)  Transformation  and  loading  of  N  and  P  species  in 
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the  overland  flow;  and  3)  Simulation  of  in-stream  water  quality  based  on 
the  QUAL2E  algorithms. 

Given  improvements  in  our  understanding  of  channel  kinetics  since  the 
release  of  QUAL2E,  this  development  effort  has  focused  on  improving  the 
NSM  channel  kinetics  to  take  advantage  of  state-of-the-art  process 
descriptions.  Stream  water  quality  can  be  dealt  with  at  different  levels  of 
detail.  The  in-stream  nutrient  kinetic  algorithms  in  NSM  Version  1.2  are 
adapted  from  the  QUAL2E,  QUAL2K,  CE-QUAL-RIVi  (RIVi),  CE-QUAL- 
W2  (W2),  and  CE-QUAL-ICM  (ICM)  models.  Both  QUAL2  and  RIVi  are 
one-dimensional  (lD)  comprehensive  stream  water  quality  models.  W2  is 
a  water  quality  and  hydrodynamic  model  in  2D  (longitudinal-vertical)  for 
rivers,  estuaries,  lakes,  reservoirs,  and  river  basin  systems  (Cole  and  Wells 
2003).  W2  models  basic  eutrophication  processes  such  as  temperature, 
nutrients,  algae,  dissolved  oxygen,  organic  matter  and  sediment  relation¬ 
ships.  ICM  is  a  finite  volume  eutrophication  model  that  can  be  used  to 
simulate  1-,  2-,  3-dimensional  water  quality  variables  including  multiple 
forms  of  algae,  carbon,  nitrogen,  phosphorus,  and  silica;  and  dissolved 
oxygen  (Cerco  and  Cole  1995).  The  QUAL2  and  RIVi  models  lack  a  clear 
operational  definition  of  water  quality  parameters  within  the  model.  For 
example,  it  is  known  that  there  are  many  forms  of  organic  nitrogen 
present  in  natural  waters.  Both  models  combine  all  of  them  under  ‘organic 
N’  and  do  not  specify  further  whether  they  are  TON,  Kjeldahl  N,  partic¬ 
ulate,  dissolved,  or  other.  In  contrast,  the  ICM  model  contains  a  precise 
specification  of  the  nitrogen,  phosphorus,  and  carbon  variables.  Since  ICM 
contains  more  processes,  it  also  contains  more  parameters.  In  principle, 
the  larger  number  of  parameters  in  models  entails  more  difficulties  during 
model  calibration. 

This  report  describes  the  in-stream  water  quality  module  within  NSM 
version  1.2.  The  in-stream  water  quality  module  includes  the  major 
processes  for  nitrogen,  phosphorus,  carbon,  and  oxygen  cycling  in  a 
stream.  Kinetic  process  equations  are  presented  for  each  state  variable  for 
water  quality  modeling  in  streams.  These  equations  were  taken  largely 
from  the  QUAL2K,  ICM,  etc.  and  summarize  past  development  efforts.  As 
research  continues,  it  is  anticipated  that  improved  process  descriptions 
will  be  developed  and  will  be  integrated  into  the  NSM. 
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2  Water  Quality  Transport  Mass  Balance 
Equations 


Water  quality  in  estuaries,  lakes,  and  reservoirs  depends  upon  the  quantity 
and  quality  of  inflows  from  the  upstream  watershed,  which  is  usually  the 
most  significant  source  of  pollutants.  Streamflow  is  therefore  of  primary 
importance  to  the  health  of  receiving  waters  and  its  quantification  is 
crucial  to  the  ability  to  manage  those  systems  in  an  environmentally 
healthy  manner. 

For  stream  water  quality  studies,  it  is  assumed  that  longitudinal  and 
temporal  changes  (lD  transport)  are  applicable.  Water  quality  is  affected 
in  streams  due  to  physical  transport  and  exchange  processes  and  bio¬ 
logical,  chemical,  and  biochemical  kinetic  processes  along  with  changes 
due  to  benthic  sediments.  Transport  not  only  includes  advection  but 
diffusion/dispersion  of  water  quality  constituents  as  well.  Derivation  of 
the  in-stream  transport  equations  can  be  found  in  Appendix  A. 


The  following  lD  transport  equation  is  often  solved  for  dissolved  state 
variables: 


dc,._ 

QdCj  l  d 

dt 

A  dx  A  dx 

dc. 

AD,  ' 


dx 


XR,* 


(1) 


where: 

Cj  =  concentration  of  dissolved  constituent  j  in  stream  [M/L3] 

Q  =  stream  flow  discharge  [L3/T] 

A  =  cross-sectional  area  of  the  stream  [L2] 

Dx  =  longitudinal  dispersion  coefficient  of  constituent  j  [L2/T] 
qi  =  volumetric  lateral  inflow  [L2/T] 
hch  =  hydraulic  depth  of  the  stream  [L] 

Cj,i  =  inflow  concentration  of  dissolved  constituent  j  [M/L3] 

Sd  =  source/sink  term  due  to  transfer  from  sediment  pore  water 
[M/L2/T] 

Csj  =  concentration  of  particulate  constituent)  in  stream  [M/L.3] 
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Csj,i  =  inflow  concentration  of  particulate  constituent  j  [M/L3] 

Ss  =  sink  term  due  to  settling  from  water  column  to  bed  sediment 
[M/L2/T] 

Sr  =  source  term  due  to  resuspension  of  bed  sediment  [M/L2/T] 
Rj,k  =  source/sink  term  representing  a  mass  rate  of  change  of 
constituent  j  due  to  kinetic  reaction  k  [M/Ls/T] 

S  j=  source/sink  term  representing  external  gains  and  losses  of 
constituent  j  [M/L3/T]. 


The  transport  of  particulate  pollutants  includes  the  settling,  resuspension, 
and  sedimentation  of  solids.  Water  quality  constituents  sorbed  onto  solid 
particles  are  transported  between  the  water  column  and  the  sediment  bed. 
Due  to  its  solid  nature,  there  is  no  transfer  from  sediment  pore  water.  The 
lD  transport  equation  for  particulate  state  variables  is: 


dCsj 

Q 

dCsj 

1 

d  ' 

dt 

A 

dx 

^  A 

dx 

dCsi 

AD .  9 


dx 


q‘  (.c„-cj~is,-sr)+ZRj*  +s, 


"  ^sj,l  ~  s  j 


ch 


(2) 


Transport  Equations  l  and  2  compute  the  transport  of  constituents  with 
their  internal  and  external  source/sink  terms.  External  terms  are  due  to 
boundary  sources/sinks  and  transfer  mechanisms.  Internal  source/sink 
terms  represent  kinetic  reaction  rates.  The  division  of  terms  allows  kinetic 
sources/sinks  to  be  updated  at  different  frequencies  than  the  transport 
terms  -  consistent  with  the  coarser  time  steps  associated  with  biological 
and  chemical  processes. 

In  the  following  mass  balance  equations,  the  transport  operator  is  defined 
as: 


L(C  )  =  ——^A 
J  A  dx 


1  d 


A  dx 


dCA 

AD  J 


dx 


+  ^-(C,,-CJ)  +  S, 


(3) 


The  exchange  of  dissolved  mass  between  the  stream  channel  and  the 
benthic  sediment  is  modeled  as  a  first-order  mass  transfer  process. 

Settling  and  resuspension  of  each  component  of  the  nutrient  pool  is 
accomplished  using  the  particle  sedimentation  and  resuspension  equation. 
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The  NSM  in-stream  water  quality  module  simulates  nutrient  (nitrogen, 
phosphorus,  and  carbon)  cycling,  dissolved  oxygen  (DO)  dynamics, 
phytoplankton  production  and  loss.  The  mass  balance  equations  for  each 
NSM  state  variable  are  presented  in  the  following  sections.  In  order  to 
provide  a  complete  description  of  the  dominant  pools  and  fluxes  in  the 
water  column,  a  single  benthic  sediment  layer  is  included  to  maintain 
mass  balance.  The  bed  sediment  plays  an  important  role  in  the  transport 
and  fate  of  water  quality  constituents.  Sediment-sorbed  pollutants  may  be 
buried  in  the  bed  by  deposition  and  sedimentation,  or  they  may  be 
released  to  the  water  column  by  scour.  Within  NSM  the  governing 
equations  are  greatly  simplified,  as  the  processes  of  advection,  dispersion, 
and  lateral  inflow  are  not  modeled  for  the  benthic  sediment.  The  concen¬ 
trations  of  the  dissolved  and  particulate  state  variables  in  a  stationary 
upper  bed  are  given  by: 


dt 


E 


(4) 


My 

dt 


(5) 


where: 

C/2  =  concentration  of  dissolved  constituent  j  in  bulk  volume  of 
benthic  sediment  [M/L3] 

Csj2  =  concentration  of  particulate  constituent)  in  benthic  sediment 
[M/L3] 

hsed  =  depth  of  upper  bed  layer  [L]. 
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3  Kinetic  Reaction  Coefficients 

A  principal  component  in  the  mass  balance  transport  equations  written  for 
nutrients  is  the  nutrient  uptake  kinetics  associated  with  algae  (phyto¬ 
plankton)  growth.  To  specify  the  nutrient  uptake  kinetics  with  this  growth 
it  is  necessary  to  specify  the  population  stoichiometry  in  units  of  nutrient 
uptake  per  mass  of  population  synthesized.  Using  carbon  as  the  unit  of 
population  biomass,  the  relevant  ratios  are  the  mass  of  nitrogen  and 
phosphorus  per  unit  mass  of  carbon.  A  selection  of  these  ratios  is 
described  as  follows. 

Stoichiometric  coefficients  of  organic  matter 

The  mass  ratios  of  carbon  to  nitrogen  to  phosphorus  to  algae  are  required 
by  the  following  mass  balance  equations.  The  stoichiometric  coefficients  of 
algae  are  calculated  with  the  technique  provided  in  Chapra  (1997),  based 
on  a  Redfield  composition  of  all  organic  particles  (Table  1).  These  recom¬ 
mended  values  are  then  combined  to  determine  stoichiometric  coefficients 
as  follows. 


r  =  7.2 


mgN 

mgA 


r  i.oESN 

mgA 


'•„  =  0.04 

mgA 


r-2'67f 


r  =  r 

cca  ca 


gC 


mgA 


x- 


moleC  m3 


-x- 


12  gC  1000  L 


r  = 

cco 


1 

gC 

moleC 

roc 

[g°J 

12  gC 

-X- 


nr 


(6) 

(7) 

(8) 

(9) 

(10) 


(11) 
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Table  1.  Recommended  values  for  stoichiometry 


Variable 

Symbol 

Recommended  values 

Carbon 

C 

40  gC 

Nitrogen 

N 

7.2  gN 

Phosphorus 

P 

lgP 

Chlorophyll  a 

A 

IgA 

Stoichiometric  coefficients  for  oxygen  generation  and  consumption 

NSM  requires  that  the  rates  of  oxygen  generation  and  consumption  be 
prescribed.  The  stoichiometric  coefficients  for  oxygen  generation  and 
consumption  are  based  upon  a  typical  chemical  reaction  for  plant  photo¬ 
synthesis  and  respiration  (Chapra  1997).  Photosynthesis  is  one  of  the  most 
important  biological  activities  in  stream  systems.  Many  water  quality  para¬ 
meters  such  as  dissolved  oxygen,  carbon  dioxide,  pH  cycles,  and  nutrients 
are  regulated  by  the  photosynthetic  reaction  in  phytoplankton.  Simply 
stated,  photosynthesis  is  the  process  by  which  phytoplankton  uses  sunlight 
to  convert  carbon  dioxide  into  a  food  source  and  to  release  oxygen  as  a  by¬ 
product.  Note  that  the  equation  for  respiration  (R)  is  the  opposite  of 
photosynthesis  (P)  indicating  the  products  of  photosynthesis  become 
reactants  in  respiration  and  vice-versa.  Because  it  requires  light,  photo¬ 
synthesis  occurs  only  during  daylight  hours.  Respiration,  on  the  other 
hand,  occurs  24  hours  a  day. 

Ammonium  as  substrate 

106CO,  +{l6NH+  +HPOl}  +  108U2O^{Clo6H263O110N16P1}  +  107O2  +  14H+ 

R 

nutrients  algal  protoplasm 

Nitrate  as  substrate 

106CO,  +  {l6AT03  + //P042  }  +  122H2O+i8H+^>{C106//263O110AT16P1}  + 13802 

R 

Oxygen  is  produced  during  photosynthesis  and  consumed  during  respi¬ 
ration.  If  ammonium  is  the  substrate,  the  following  ratio  can  be  used  to 
determine  the  grams  of  oxygen  generated  for  each  gram  of  plant  matter 
that  is  produced  through  photosynthesis. 
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107  mole02(32  g02/mole02)  _  2  6Qg02 
106  moleC(12  gC/moleC)  gC 


(12) 


If  nitrate  is  the  substrate,  the  following  ratio  applies: 

r  138  mole02  (32  g02/mole02 )  3  ^  ?  g02 

106  moleC(12  gC/moleC)  gC 

Nitrification  and  denitrification 

NH+  +  202  >  NO:  +  H„0  +  2H+ 

4  2  3  2 

5CH20  +4NO3  +4H+  denitrification  >  +  ^  +^0 


Nitrification  is  dependent  upon  a  suitable  supply  of  DO.  The  ratio  of 
oxygen  to  nitrogen  in  the  nitrification  process  is  determined  by: 

r  _  2  mole02(32  g02/mole02)  _  4  5?g02 
1  moleN  (14  gN/moleN)  gN 


Stoichiometric  coefficient  for  DOC  utilization 

DOC  is  utilized  during  denitrification.  The  following  ratio  ( rcndn )  is  used  to 
determine  the  organic  carbon  lost  per  nitrate  nitrogen  that  is  denitrified. 

r  5  moleCx  12  gC/moleC  x  igN  1Q711  §c  (15 
cndn  4  moleN  x  14  gN/moleN  1000  mgN  mgN 


Temperature  effects  on  kinetic  reaction  rates 

All  rates  depend  exponentially  on  temperature.  The  temperature  effect  for 
all  first-order  reactions  used  in  the  model  is  represented  by: 

k(T)  =  k(20)0T~20  (16) 


where: 

k{T)  =  reaction  rate  at  temperature  T°C  [T1] 
k{ 20)  =  reaction  rate  at  temperature  20°C  [T1] 

0  =  temperature  coefficient  for  the  reaction. 
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4  Water  Quality  Kinetic  Reaction  Equations 

Well-developed  models  and  tools  are  available  to  address  the  physical 
transport  components  of  equations  in  water  quality  modeling.  Particularly, 
lD,  2D,  and,  increasingly,  3D  hydrodynamic  models  are  available  to 
determine  the  velocity  field  and  are  becoming  more  practical  with 
advancements  in  computer  technology.  Therefore,  future  water  quality 
modeling  developments  will  refine  the  description  of  kinetic  processes  in 
Equations  1  and  2.  Kinetic  processes  describe  changes  in  the  constituent 
concentrations  that  are  due  to  biological,  chemical,  biochemical,  and 
physical  processes.  The  historical  development  of  oxygen,  nitrogen, 
phosphorus,  and  carbon  water  quality  models  shows  step-by-step 
extensions  and  increasing  complexity. 

This  chapter  describes  governing  equations  for  the  current  in-stream 
water  quality  module,  which  includes  the  major  processes  for  nitrogen, 
phosphorus,  carbon,  and  oxygen  cycling.  The  water  quality  processes  were 
simplified  on  the  basis  of  selecting  the  dominant  biochemical  processes 
under  consideration.  Only  process  descriptions  applicable  to  small 
streams  are  being  incorporated  into  NSM.  Hydrolysis  of  particulate 
organic  nutrients  into  dissolved  organic  form  is  modeled  as  a  first-order, 
temperature-dependent  process.  Mineralization  of  the  organic  nutrient 
pools  back  to  inorganic  nutrients  (i.e.  dissolved  inorganic  ammonium, 
phosphorus,  carbon)  is  also  modeled  as  a  first-order,  temperature- 
dependent  process.  The  terms  for  mortality  and  excretion  are  modeled  in 
the  usual  manner  as  first  order  losses  with  an  identical  temperature 
dependence  for  growth.  Nitrification  and  denitrification  influence  the 
dissolved  inorganic  nitrogen  fraction,  which  includes  NHA  and  N03  as 

state  variables.  These  processes  are  represented  as  first  order,  oxygen- 
dependent  reactions.  The  NSM  in-stream  water  quality  module 
incorporates  12  state  variables  (Table  2). 
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Table  2.  In-stream  water  quality  state  variables 


Variable 

Symbol 

Notes 

Nitrate  Nitrogen 

NO - 

Ammonium  Nitrogen 

NH+ 

Dissolved  Organic  Nitrogen 

DON 

Particulate  Organic  Nitrogen 

PON 

Dissolved  Inorganic  Phosphorus 

DIP 

Dissolved  Organic  Phosphorus 

DOP 

Particulate  Organic  Phosphorus 

POP 

Dissolved  Inorganic  Carbon 

DIC 

Dissolved  Organic  Carbon 

DOC 

Particulate  Organic  Carbon 

POC 

Dissolved  Oxygen 

DO 

Algae  1 

Api 

Phytoplankton 

Algae  2 

Ap2 

Phytoplankton 

Algae  3 

Abi 

Bottom  algae 

Algae  4 

Ab2 

Bottom  algae 

These  model  state  variables  can  be  used  to  compute  the  following 
composite  variables: 

•  Total  Nitrogen  (TN) 

•  Total  Kjeldahl  Nitrogen  (TKN) 

•  Total  Phosphorus  (TP) 

•  Total  Organic  Carbon  (TOC) 

•  Ultimate  Carbonaceous  BOD  (CBOD) 

TN  (TP)  is  a  measure  of  all  forms  of  dissolved  and  particulate  nitrogen 
(phosphorus)  present  in  waters.  TOC  is  a  measure  of  all  the  various  forms 
of  organic  C  found  in  waters.  Ultimate  CBOD  reflects  oxidation  of  both 
dissolved  and  particulate  organic  carbon. 

TN  =  N03  +  NH;  +  PNH4+  +  DON  +  PON 


TKN  =  NHa  +  PNH+a  +  DON  +  PON 
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TP  =  DIP  +  PIP  +  DOP  +  POP 
TOC  =  DOC  +  POC 
CBOD  =  rocDOC  +  rocPOC 

The  transport  of  NSM  state  variables  with  their  kinetic  reaction  rates 
expressed  in  source  and  sink  terms  is  described  in  the  following  sections. 
Additional  state  variables  such  as  benthic  algae,  zooplankton,  and 
temperature  will  be  added  in  the  future  versions  of  NSM. 

Nitrogen 

Nitrogen  is  transported  into  stream  systems  from  airborne,  land  surface, 
underground,  and  in  situ  sources.  Much  of  the  nitrogen  entering  streams 
is  transported  with  eroded  sediments  and  solid  organic  matter  (SOM),  and 
dissolved  form  in  surface  runoff.  The  nitrogen  forms  commonly  found  in 
river  water  are:  N03 ,  N02 ,  NH4 ,  DON,  and  PON  (Meybeck  1982).  These 

forms  are  reactive  in  the  framework  of  the  N  cycle.  The  nitrogen  cycle 
includes  the  additional  processes  of  denitrification,  nitrification,  and  N2 
fixation  that  are  not  in  the  carbon  and  phosphorus  cycles.  N02 

concentrations  are  not  tracked  by  the  model  because  the  amount  of 
N02  present  in  the  stream  is  usually  very  small  relative  to  N03  .  The  in- 

stream  nitrogen  cycle  simulated  in  NSM  is  illustrated  in  Figure  1. 

The  major  processes  involved  in  nitrogen  transformations  in  the 
sediments  and  water  column  are: 

•  Mineralization  of  DON  to  NH4 

•  Nitrification  of  NH4  to  N03 

•  Denitrification  of  N03 

•  Biological  uptake  of  NH4  and  N03  by  phytoplankton  and  algae 

•  Biological  respiration  of  NHA  by  phytoplankton  and  algae 

•  Dissolved  sediment  fluxes  of  NH4  ,  N03 ,  and  DON 

•  Adsorption/desorption  of  NH4  onto  suspended  sediments 

•  Hydrolysis  of  PON  to  DON 

•  Biological  mortality  and  excretion  into  the  DON  and  PON  pools 

•  Settling/resuspension  of  PON  and  sediment  attached  inorganic 
nitrogen. 
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Figure  1.  Simplified  schematic  representation  of  the  in-stream  nitrogen  cycling  processes 


Four  nitrogen  variables  are  modeled:  PON,  DON,  NH+ ,  and  AT03 . 

Ammonium  nitrogen  is  divided  into  particulate  and  dissolved 
concentrations  by  spatially  variable  dissolved  fractions  reflecting  sorption. 

The  nitrogen  equations  in  the  water  column  and  benthic  sediment  are 
summarized  as  follows: 


Water  column 


B  PON  1 

— =  UPONa)  +  —(vsPONa  -  v,.PON 
at  h. 


sed  ) 


Lch 


+  rnaKPAP  +  J^QonKiA,  ~  KnPONch 


Lch 


(17) 


=  L(DON,h)  +  ^(DON,Jcp- DON A 
9t  hch 

+  hnPONch-FoxmnkmnDONch 


(18) 


BTAfl-f+  k 

- ^ -  =  UTNH;ct )  +  j-i-( fd2TNH;xJ<p  -  fdTNH )  +  FoxmnkmnDON ch 


ch 


h 


ch 


+ 

4  sed 


n=l 


71=1 


(19) 


rnaKpAp  ~  rnaPapBpAp  ~  PabPbNAb  ~  FoxnaknitfdTNH . 


4  ch 


ch 
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9N03ch 

dt 


UNO-3A ) + i(J -  no;a  ) + F0Jc^fdTNH;,h 

nch 

- (1  -  F )-  D0Cntn  k„NO;A  - r„„ (1  -Pv)ntAp 

Kscdn+  DULsed 

~  ^  —  Pab  ^bN^-b 


(20) 


Benthic  sediment 


dPONsed 

dt 


( v,PON , 


O'PON^-k^PON 


sed 


(21) 


dDON 


sed 


dt 


h 


DONch-DONsed 


M- 


khnPONsed-F, 


kDONspd  (22) 


oxmn  mn 


sed 


4  -  '  [f,TNH*,A  -  fd2TNH dxd /<p 


dt 


h 


sed 


Faxmnb . DONmd  -  FmknJd2TNH4 

-'tfr2dvrdTNH 


h 


sed 


+ 

4  sed 


(23) 


dNO~sed 

~  K 

dt 

Ked 

N0-3di-^ML 


(A  F  nxdn  ) 


(P 

DOC 


FOXnaKitfd2TNH;sed 


sed 


F-scdn  +  F>OCsed 


KnitNO,  sed 


(24) 


where: 


PONch 

Vs 

Vr 

kdp 

kdb 

AP 

Ab 

BA 

khn 


concentration  of  in-stream  PON  [M/Ls] 
particle  settling  rate  [L/T] 
particle  resuspension  rate  [L/T] 

temperature-dependent  phytoplankton  death  rate  [T1] 
temperature-dependent  bottom  algae  death  rate  [T1] 
stream  phytoplankton  concentration  [M/L3] 
stream  bottom  algae  concentration  [M/L2] 
stream  bottom  surface  area  [L2] 

temperature-dependent  PON  hydrolysis  rate  coefficient  [T1] 
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DONch  =  concentration  of  in-stream  DON  [M/L3] 

kmn  =  temperature-dependent  DON  mineralization  rate  coefficient 
[T1] 

Foxmn  =  DON  mineralization  attenuation  due  to  low  oxygen 
TNH^ch  =  total  concentration  of  in-stream  NH'A  [M/L3] 

krp  =  temperature-dependent  phytoplankton  respiration  rate  [T1] 
Foxna  =  nitrification  attenuation  due  to  low  oxygen 
knit  =  temperature-dependent  NHf  nitrification  rate  coefficient  [T1] 

Pap  =  preference  coefficient  of  phytoplankton  for  NHA 
Pab  =  preference  coefficient  of  bottom  algae  for  NHA 

pp  =  phytoplankton  photosynthesis  rate  [T1] 

\ib  =  bottom  algae  photosynthesis  rate  [T1] 

N03ch  =  concentration  of  in-stream  N03  [M/L.3] 

Kscdn  =  DOC  half-saturation  constant  for  denitrification  [gC/m3] 
[M/L3] 

kdnit  =  temperature-dependent  jV03  denitrification  rate  coefficient 
[T1] 

Foxdn  =  effect  of  low  oxygen  on  denitrification 
fPn  =  fraction  of  the  total  concentration  in  particulate  form 

associated  with  particle  “n” 

fd  =  fraction  of  the  total  concentration  in  dissolved  form 
PONsed  =  concentration  of  bed  sediment  PON  [M/L3] 

DONsed  =  concentration  of  bed  sediment  DON  [M/L3] 

(j)  =  sediment  porosity 
Ke  =  mass  transfer  coefficient  [M/T] 

TNH^sed  =  total  concentration  of  bed  sediment  NHA  [M/L.3] 

N03sed  =  concentration  of  bed  sediment  AT03  [M/L3], 

Phosphorus 

Phosphorus  is  an  important  element  because  it  is  usually  in  short  supply 
relative  to  other  macronutrients  (Correll  1998;  Carpenter  et  al.  1998).  In 
contrast  to  nitrogen  and  carbon,  there  is  no  gaseous  atmospheric  source. 
Phosphorus  enters  rivers  primarily  as  particulate  matter  and  secondarily 
as  DIP  also  known  as  “ortho-P”  (H3P04  and  conjugate  base  forms).  The 
phosphorus  forms  simulated  in  NSM  are  DIP,  DOP,  and  POP.  The  in- 
stream  phosphorus  cycle  simulated  in  NSM  is  illustrated  in  Figure  2. 
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The  major  processes  involved  in  phosphorus  transformations  in  the 
sediments  and  water  column  are: 

•  Mineralization  of  DOP  to  DIP 

•  Biological  uptake  and  respiration  of  DIP  by  phytoplankton  and  algae 

•  Dissolved  sediment  fluxes  of  DIP  and  DOP 

•  Adsorption/desorption  of  DIP  onto  suspended  sediments 

•  Hydrolysis  of  POP  to  DOP 

•  Biological  mortality  and  excretion  into  the  DOP  and  POP  pools 

•  Settling/resuspension  of  POP  and  PIP. 


nminim 


Figure  2.  Simplified  schematic  representation  of  the  in-stream  phosphorus  cycling  processes 

Three  phosphorus  variables  are  modeled:  POP,  DOP,  and  inorganic 
phosphorus.  Inorganic  phosphorus  is  divided  into  particulate  and 
dissolved  concentrations  by  spatially  variable  dissolved  fractions. 
Dissolved  or  available  inorganic  phosphorus  (DIP)  interacts  with 
particulate  inorganic  phosphorus  (PIP)  via  a  sorption-desorption 
mechanism. 

The  phosphorus  equations  in  the  water  column  and  benthic  sediment  are 
summarized  as  follows: 

Water  column 


dPOPch 

dt 


--U.POPA)  +  ^-[v,POPA 

nch 

1 

+  rpakdpAp+^-<l0PkdbAb 

•Lnh 


»r POP „) 
KnPOPA 


(25) 
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dDQPch 

dt 


~~  UDOPct )  +  -r~{DOPsedl<p -  DOPd, ) 

nch 

+  khpPOPrfi  ~  Poxmp^mpPPPch 


(26) 


9TIPch 

dt 


L(TIPch )  +  ^{fd2TIPsed /cp  -  fdTIPch ) 


Lch 


ch 


N  N 


+  Kmp KpDOPj,  +  rpXpAp  ~  ,-4,.  -  -j^VtrA 


ch 


(27) 


Benthic  sediment 


dPOPsed 

dt 


{ vPOP 


ch 


vrpOPsed)-khpPOPsed 


(28) 


ODOP 


sed 


dt 


h 


(DOPct  -  DOP„Jcp)  +  khpPOP£ed  -  FopmpkmpDOPxi  (29) 


sed 


dTIP. 


sed 


dt 


h 


sed 

-  FoxmpKpDOPsed 


sed 


N 


N 


—(fJipA-fi2TipxdM+-^-\Y,fppP,TIppi,  -E  fp2pPrTIP^ 


\n= 1 


n= 1 


(30) 


where: 


POPch 
DOPch 
TIPch 
POPsed 
DOPsed 
TIP sed 
khp 
kmp 

Foxmp 

fd‘2 

fP2n 


concentration  of  in-stream  POP  [M/Ls] 

concentration  of  in-stream  DOP  [M/Ls] 

total  concentration  of  in-stream  inorganic  P  [M/L.3] 

concentration  of  bed  sediment  POP  [M/L.3] 

concentration  of  bed  sediment  DOP  [M/L.3] 

total  concentration  of  bed  sediment  inorganic  P  [M/L3] 

temperature-dependent  POP  hydrolysis  rate  coefficient  [T1] 

temperature-dependent  DOP  mineralization  rate  coefficient 

[T1] 

DOP  mineralization  attenuation  due  to  low  oxygen 
fraction  of  total  concentration  in  the  dissolved  phase 
fraction  of  total  concentration  in  the  particulate  phase 
associated  with  size  fraction. 
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Carbon 

The  carbon  constituents  taken  into  account  are:  particulate  organic  carbon 
(POC),  dissolved/colloidal  organic  carbon  (DOC),  and  dissolved  inorganic 
carbon  (DIC).  DOC  in  streams  is  usually  comprised  of  minor  amounts  of 
biodegradable  plant,  phytoplankton  and  bacterial  residues,  which  are 
being  rapidly  recycled,  and  major  amounts  of  biological  refractory  residues 
comprised  of  substituted  aromatic  structures,  and  branched  and  partially 
oxidized  cyclic  aliphatic  structures.  The  DIC  actually  consists  of  several 
carbonate  constituents:  carbon  dioxide  (C02),  bicarbonate  (HC03),  and 
carbonate  (C032).  The  value  of  water  pH  is  controlled  by  the  carbonate 
system.  Within  the  stream  POC  is  not  inert  but  readily  transformed  into 
DIC  either  directly  or  via  DOC  as  an  intermediate  step.  In-stream  DOC  is 
produced  by  incomplete  hydrolysis  of  POC  or  by  desorption  from  mineral 
surfaces.  DOC  contributes  to  the  turbidity  of  surface  waters,  provides  an 
energy  source  for  organisms,  and  affects  the  fate  and  bioavailability  of 
contaminants.  Hydrophobic  chemicals  and  some  metals  have  high  binding 
affinities  for  DOC.  The  dynamics  of  these  toxics  is  intimately  connected 
with  the  generation,  transport,  and  fate  of  organic  carbon.  Nitrogen  and 
phosphorus  are  contained  as  nutrients  within  DOC.  The  DIC  balance 
additionally  includes  atmospheric  fluxes  of  C02  that  are  based  on  the 
difference  between  the  atmospheric  and  water  column  values  of  pC02 
(vapor  pressure  of  carbon  dioxide).  To  estimate  the  C02  fraction  of  the 
DIC  pool,  the  carbonate  buffer  system,  alkalinity,  and  pH  that  govern  the 
subsequent  partitioning  of  DIC  between  pC02,  bicarbonate,  and  carbonate 
ions  need  to  be  modeled.  The  gaseous  and  aqueous  phase  C02  values  are 
related  by  Henrys  Law  for  calculation  of  pC02.  The  in-stream  carbon  cycle 
simulated  in  NSM  is  illustrated  in  Figure  3. 

The  major  processes  involved  in  carbon  transformations  in  the  sediments 
and  water  column  are: 

•  Atmospheric  fluxes  of  DIC 

•  Mineralization  of  DOC  to  DIC 

•  Biological  uptake  and  respiration  of  DIC  by  phytoplankton  and  algae 

•  Dissolved  sediment  fluxes  of  DIC  and  DOC 

•  Hydrolysis  of  POC  to  DOC 

•  Biological  mortality  and  excretion  into  the  DOC  and  POC  pools 

•  Settling/resuspension  of  POC. 
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Figure  3.  Simplified  schematic  representation  of  the  in-stream  carbon  cycling  processes 


Three  carbon  variables  are  modeled:  POC,  DOC,  and  DIC.  The  carbon 
equations  in  the  water  column  and  benthic  sediment  are  summarized  as 
follows: 

Water  column 
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nch 

+  rccoFoXmckmcD°Cch  +  PcaKPAp  -Vcca\XpAp 

1 

+  l-(.rccaKbAb-rccaVbAb^ 
nch 


(33) 
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POCch  =  concentration  of  in-stream  POC  [M/L3] 

DOCch  =  concentration  of  in-stream  DOC  [M/Ls] 

DICch  =  concentration  of  in-stream  DIC  (mole/L)  [M/L.3] 

POCsed  =  concentration  of  bed  sediment  POC  [M/L3] 

DOCsed  =  concentration  of  bed  sediment  DOC  [M/Ls] 

DICsed  =  concentration  of  bed  sediment  DIC  (mole/L)  [M/L3] 

khc  =  temperature-dependent  POC  hydrolysis  rate  coefficient  [T1] 
Foxmc  =  DOC  mineralization  attenuation  due  to  low  oxygen 
kmc  =  temperature-dependent  DOC  mineralization  rate  [T1] 
kac  =  0.923  ka  =  temperature-dependent  C02  reaeration  coefficient 
[T1] 

kH  =  Henry's  constant  [mole/(L  •  atm)] 

Pco2  =  Partial  pressure  of  carbon  dioxide  in  the  atmosphere  [atm] 


ao  =  fraction  of  total  inorganic  carbon  in  carbon  dioxide. 


Algae  Group 


Algae  are  commonly  included  as  state  variables  in  water  quality  models 
because  they  impact  DO  and  material  cycling  in  water  bodies  and  because 
excessive  algae  populations  are  of  environmental  concern.  Phytoplankton 
and  bottom  algal  communities  are  characterized  as  part  of  a  physical, 
chemical,  and  biological  assessment  of  water  quality.  Algae  that  are 
suspended  or  actively  swimming  in  the  water  column  are  referred  to  as 
phytoplankton.  The  term  "periphyton"  is  sometimes  used  to  refer  to 
benthic  algae,  while  others  use  the  term  to  refer  to  the  entire  attached 
community  of  microorganisms,  including  algae,  bacteria,  fungi,  and 
protozoa.  To  avoid  confusion  within  this  section,  “bottom  algae”  are  used 
to  designate  the  algal  community  attached  to  bottom  rocks  and  stable  sand 
surfaces.  In  many  shallow  streams  and  rivers,  it  is  the  attached  algae 
(benthic  algae,  or  periphyton,  attached  to  submerged  substrates)  that  are 
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often  of  greater  importance.  These  attached  plants  affect  water  quality  in 
various  ways,  and  their  impact  must  often  be  considered  in  order  to 
properly  evaluate  riverine  water  quality  conditions. 

NSM  simulates  algae  that  float  in  water  bodies  for  phytoplankton.  Algal 
involvement  in  the  nutrient  cycles  is  depicted  schematically  in  the 
previous  sections.  Algae  reduce  the  concentration  of  nutrients  in  the  water 
by  consuming  phosphates,  nitrate,  and  ammonium.  Through  assimilation 
these  nutrients  are  transformed  into  organic  materials  which  serve  as  a 
food  source.  A  portion  of  the  organic  matter  that  is  not  used  for  food 
decomposes,  which  further  affects  the  oxygen  and  nutrient  levels  in  the 
water  bodies.  The  activities  of  algae  in  turn  are  affected  by  the  physical 
environment.  Figure  4  provides  an  overview  of  algal  kinetic  processes 
simulated  in  NSM. 
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Figure  4.  Simplified  schematic  representation  of  the  in-stream  algae  kinetic  processes 

NSM  simulates  the  growth  and  accumulation  of  two  forms  of  bottom 
algae.  The  bottom  algae  algorithm  includes  algae  biomass  as  a  state 
variable  and  three  multi -variate  rate  vectors:  primary  production,  algal 
respiration,  and  the  removal  processes  that  result  in  export  of  biomass 
from  the  system. 

Four  algae  groups  are  configurable  within  NSM.  Algal  biomass  can  be 
simulated  either  in  units  of  chlorophyll-a  (pg  Chla  L_1)  or  carbon 
(mg  C  L-1).  Chlorophyll-a  is  most  commonly  available  as  an  estimate  of 
algal  biomass.  Algal  biomass  can  be  represented  in  the  model  either  for  the 
entire  algae  assemblage  or  as  the  individual  contributions  by  certain 
groups,  e.g.,  phytoplankton,  benthos,  and  periphyton. 
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Phytoplankton 

The  phytoplankton  biomass  increases  due  to  photosynthesis  and  is  lost  via 
respiration,  death,  and  settling.  Through  the  process  of  photosynthesis, 
phytoplankton  consume  carbon  dioxide  and  release  oxygen  back  into  the 
water  body.  At  the  same  time  algal  respiration  consumes  oxygen  and 
releases  carbon  dioxide. 

NSM  allows  for  the  specification  of  two  phytoplankton  groups.  Each  is 
governed  by  the  same  mass  balance  equation.  Distinctions  between  the 
groups  are  represented  by  using  different  parameter  values.  The  algo¬ 
rithms  for  predicting  variations  in  phytoplankton  concentrations  are  based 
upon  routines  included  in  the  QUAL2K  model  (Chapra  et  al.  2005).  The 
kinetic  equation  for  in-stream  phytoplankton  can  be  expressed  as: 


dA  1 

=  -fr-{vsAp  -VrApsed)  +  HpAp-kn>Ap -KPAP  (37) 

Without  interactions  with  the  settled  algae  flux,  the  kinetic  equation  for  in- 
stream  phytoplankton  can  be  reduced  to: 

dA  1 

-^  =  —j^vsAP+^PAp-KpAP  -Kpap  (38) 

Photosynthesis  rate 

The  photosynthetic  process  is  important  in  producing  oxygen.  Nutrients 
tend  to  limit  photosynthesis  with  phosphorus  commonly  the  limiting 
nutrient  due  to  its  low  mobility. 

Phytoplankton  photosynthesis  is  a  function  of  temperature,  nutrients,  and 
light: 


A Lp  =  KpVNpVLp 


(39) 


where: 


pp  =  phytoplankton  photosynthesis  rate  [T1] 
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kgp  =  maximum  photosynthesis  rate  at  temperature  T  [T1] 

(j )np  =  phytoplankton  nutrient  attenuation  factor  (between  o  and  l) 

(j )lp  =  phytoplankton  light  attenuation  coefficient  (between  o  and  l). 


Nutrient  limitation  of  the  photosynthesis  rate  is  dependent  on  nutrient 
concentration.  Michaelis-Menten  equations  are  used  to  represent  growth 
limitation  for  inorganic  N  and  P.  The  minimum  value  is  then  used  to 
compute  the  nutrient  attenuation  factor: 


(PNp  =  mm 


nh;a+no~ 


(40) 


where: 

Ksnp  =  N  half-saturation  constant  [pg  N/L] 

Kspp  =  P  half-saturation  constant  [pg  P/L] 

KsCp  =  inorganic  C  half-saturation  constant  [mole/L]. 


The  effect  of  light  on  phytoplankton  growth  is  complicated  by  the  fact  that 
several  factors  have  to  be  integrated  to  come  up  with  the  total  effect.  It  is 
assumed  that  light  attenuation  through  the  water  follows  the  Beer- 
Lambert  law 


Iz  =  0A7Ioe~kexZ 


(41) 


where: 

Iz  =  photosynthetically  available  radiation  (PAR)  at  depth  z  below 
the  water  surface  [Ly/d] 
kex  =  light  extinction  coefficient  [L  ‘] 

I0  =  light  intensity  at  the  water  surface  [cal/cm2/d]. 

Three  models  are  used  to  characterize  the  impact  of  light  on 
photosynthesis  and  compute  the  phytoplankton  light  attenuation 
coefficient  (Table  3). 
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Table  3.  Phytoplankton  light  attenuation  equations 


No. 

Model 

Equation 

Reference 

1 

Half-saturation  model 

I  e~KxZ 

(p  —  0 

Lp  K  +1  e~KxZ 

Baly  (1935) 

2 

Smith’s  equation 

I  e~KxZ 

<Plp=  1 - 5 - y 

V*t,+(  v*-‘) 

Smith  (1936) 

3 

Steele’s  equation 

J  q~Kxz  J  q~Kxz 

<Pl,=  °  exp  1+  0 

KLp  (  KLp  J 

Steele  (1962) 

where  Klp  is  the  appropriate  phytoplankton  light  parameter  for  each  light 
model. 


Figure  5.  Three  models  used  for  phytoplankton  photosynthetic  light  dependence 


Preference  coefficient  for  ammonium  nitrogen 

Nitrogen  is  essential  to  algae  for  assimilation  of  proteins  and  enzymes. 
Although  phytoplankton  take  up  and  use  both  NHA  and  N03 ,  their 

preference  for  the  former  has  been  demonstrated  for  physiological  reasons 
(Walsh  and  Dugdale  1972).  The  preference  coefficient  of  phytoplankton 
for  NHA  as  a  nitrogen  source  is  determined  by: 


ER  DC/EL  TR-09-4 


24 


NHlhNO~ch 


ap 


hnxp  +nh:axK«d+no;J 


hnxp 


(NH:rh+NO-,hXkt^r  +  NO. 


3  ch  J 


(42) 


where: 

khnxp  =  preference  coefficient  of  phytoplankton  for  NH+  [mg  N/m3]. 


Bottom  algae 

As  with  phytoplankton,  bottom  algae  growth  is  impacted  by  temperature, 
light,  and  nutrients.  The  growth  of  bottom  algae  consumes  nutrients  and 
produces  oxygen.  Bottom  algae,  like  phytoplankton,  also  excrete  cell 
contents  and  die,  recycling  dissolved  and  particulate  organic  matter  to  the 
stream’s  carbon  and  nutrient  pools.  While  the  modeling  approaches  used 
for  phytoplankton  and  bottom  algae  are  similar,  bottom  algae  differ  from 
phytoplankton  in  a  number  of  fundamental  ways: 

•  Bottom  algae  do  not  move  with  the  water  current,  as  do  phytoplankton, 

•  Bottom  algae  typically  dwell  on  or  near  the  bottom,  so  are  not  impacted 
by  the  average  light  in  the  water  column  but  the  light  reaching  the 
bottom  (substrate). 

•  Bottom  algae  are  limited  by  the  amount  of  substrate  available  for 
growth. 

•  There  is  typically  a  maximum  density  for  attached  plants. 

Sources  and  sinks  for  bottom  algae  include  growth,  death,  and  respiration. 
Growth  is  computed  from  a  maximum  rate  that  is  then  modified  based 
upon  available  light  and  internal  nutrients.  Unlike  phytoplankton,  bottom 
light  rather  than  average  water  column  light  is  used  in  the  computation  of 
growth.  Rates  of  death  and  respiration  are  temperature  dependent.  Rates 
of  growth,  respiration,  and  death  impact  other  model  state  variables 
including  dissolved  oxygen  and  nutrients. 


The  kinetic  equation  for  in-stream  bottom  algae  can  be  expressed  as: 


dAb 


dt 


fbA>  h'oxbkrbAb  kdbAb 


(43) 
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where: 

Ab  =  stream  bottom  algal  concentration  [M/L2] 
l_i6  =  benthic  algal  photosynthesis  rate  [T1] 

Foxb  =  attenuation  due  to  low  oxygen 

krb  =  temperature-dependent  benthic  algal  respiration  rate  [T1] 
kdb  =  temperature-dependent  benthic  algal  death  rate  [T1]. 

Photosynthesis 

Bottom  algal  photosynthesis  is  a  function  of  temperature,  nutrients,  and 
light: 


A h  =  kgbVmVLb 


(44) 


where: 

kgb  =  maximum  bottom  algal  rate  at  temperature  T  [T1] 

=  bottom  algal  nutrient  attenuation  factor  (between  o  and  l) 

§>Lb  =  bottom  algal  light  attenuation  coefficient  (between  o  and  l). 

The  effect  of  nutrient  limitation  on  bottom  plant  photosynthesis  is 
modeled  in  a  different  way  than  for  the  phytoplankton.  Rather  than  being 
dependent  on  external  nutrient  concentration,  the  photosynthesis  rate  is 
dependent  on  intracellular  nutrients  using  a  formulation  originally 
developed  by  Droop  (1974) 


Tm,  =  min 


1 


QpiV  A 

Qn 


Qo P  DICch 

qp  ’KsCb  +DIC, 


(45) 


where: 

qN  and  qp  =  the  cell  quotas  of  nitrogen  [mgN  mgA”1]  and  phosphorus 
[mgP  mgA  '],  respectively, 

qo  V  and  qop  -  the  minimum  cell  quotas  of  nitrogen  [mgN  mgA”1]  and 

phosphorus  [mgP  mgA  '],  respectively,  and 
ksC  =  inorganic  carbon  half-saturation  constant  for  the  bottom 

algae  [mole/L]. 
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The  cell  quotas  represent  the  ratios  of  the  intracellular  nutrient  to  the 
bottom  plant’s  weight.  Their  mass  balance  equations  are  described  as 
follows. 

VN  =  (46a) 

Ab 

Qn  =  Ijl  (46b) 

Ab 

where: 

IN  =  intracellular  nitrogen  concentration  [mgN /m2] 

IPh  =  intracellular  phosphorus  concentration  [mgP/m2]. 

The  change  in  intracellular  nitrogen  and  phosphorus  in  bottom  algal  cells 
is  calculated  from 


dIN„ 

dt 

—  PfcAT  A>  ^oxbKb^b  Kb^b 

(47a) 

dIPt 

dt 

=  Pw>A> — ^oxbKb^b — Kb^b 

(47b) 

r\  _ 

NH;ch+N03ch  NH;chkhnxp 

(48a) 

mN  Km  +  NHtch  +  N03 ch  K N  +  (?JV  +  % 

Pw>  =  p 

DIP  kqP 

mF  Kpb  +  DIP  Kp  +  Kip  +  q0P) 

(48b) 

where: 

pmJvand  p  =  the  maximum  uptake  rates  for  nitrogen  [mgN/mgA/d] 

and  phosphorus  [mgP/mgA/d],  respectively, 
kgNb  and  kspb  =  half-saturation  constants  for  external  nitrogen  [pgN/L] 

and  phosphorus  [pgP/L],  respectively,  and 
K^n  and  =  half-saturation  constants  for  intracellular  nitrogen  [mgN 

mgA”1]  and  phosphorus  [mgP  mgA”1],  respectively. 
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Light  limitation  is  determined  by  the  amount  of  photosynthetically  active 
radiation  (PAR)  reaching  the  bottom  of  the  water  column.  This  quantity  is 
computed  with  the  Beer-Lambert  law  evaluated  at  the  bottom  of  the  river: 

IH  =  I0e~kexhck  (49) 


where: 

1(H)  is  light  intensity  at  depth  H  below  the  water  surface  [Ly/d], 

I( o)  is  light  intensity  just  below  the  water  surface  [Ly/d]. 

Three  models  are  used  to  characterize  the  impact  of  light  on  bottom  algae 
photosynthesis.  Substituting  the  above  formulation  into  these  models 
yields  the  following  formulas  for  the  bottom  algae  light  attenuation 
coefficient: 


Table  4.  Bottom  algae  light  attenuation  equations 


No. 

Model 

Equation 

1 

Half-saturation  model 

J  g—kexhch 

Lb  rr  .  -r  —lchu 

KLb  +I0e  exch 

2 

Smith’s  equation 

j  g-KxKh 

VLb=  I — 5 - j 

3 

Steele’s  equation 

T  0-kexKh  T  p-Kxhch 

m  —  °e  -1  i  ±0K 

Lb  jr  eXP  XT 

Lb  {  Lb 

where  Ku  is  the  appropriate  bottom  algae  light  parameter  for  each  light 
model. 

Losses 

Bottom  algal  biomass  decreases  due  to  respiration  and  death.  Bottom  algal 
respiration  is  represented  using  first-order,  temperature-corrected 
kinetics: 


Kb  —  ^(20)9r6 


r-20 


(50) 
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where  krb( 20)  is  the  bottom  algae  respiration  rate  constant  at  20  °C  [T  '] 
and  Orb  is  the  bottom  algae  respiration  temperature  coefficient. 

Bottom  algal  death  is  also  represented  using  first-order,  temperature- 
corrected  kinetics: 


^  =  ^(20)6,/-”  (51) 

where  kdb(20)  is  the  bottom  algae  death  rate  constant  at  20  °C  [T  ']  and  Odb 
is  the  bottom  algae  death  temperature  coefficient. 

The  preference  coefficient  of  bottom  algae  for  Nil l  as  a  nitrogen  source  is 
determined  by: 


NHl,hN03ch 


'  ab 


ft 


hnxb  +  NHtch  KKnxb  +  N03  ch  ) 

_ NH4Chkhiab _ 

( NHt  h  +  N03  h )  (kh  b  +  N03  h ) 


(52) 


where: 

khnxb  =  preference  coefficient  of  bottom  algae  for  NHA  [mg  N/m3]. 

Dissolved  Oxygen 

The  final  in-stream  cycle  that  will  be  discussed  is  the  DO  balance.  An 
adequate  DO  concentration  is  a  basic  requirement  for  a  healthy  aquatic 
ecosystem.  DO  concentration  is  generally  viewed  as  an  indicator  of  the 
overall  well-being  of  streams  and  their  associated  ecosystems.  DO 
dynamics  include  atmospheric  exchange,  the  sediment  oxygen  demand 
(SOD),  microbial  use  during  organic  matter  mineralization  and 
nitrification,  photosynthetic  oxygen  production  and  respiratory  oxygen 
consumption,  and  respiration  by  other  optional  biotic  components.  The 
major  processes  involved  in  modeling  DO  are  shown  in  Figure  6. 

The  following  processes  are  considered  in  computing  DO: 

•  Exchange  to  and  from  the  air/water  interface 

•  Utilization  of  oxygen  at  the  sediment/water  interface  (i.e.  SOD) 
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•  Utilization  of  oxygen  as  bacteria  degrade  organic  matter 

•  Utilization  of  oxygen  in  the  process  of  nitrification 

•  Photosynthetic  oxygen  production  and  respiratory  consumption  by 
phytoplankton  and  algae. 


Figure  6.  Simplified  schematic  representation  of  the  in-stream  DO  dynamic  processes 


Sources  of  DO  include  algal  photosynthesis  and  atmospheric  reaeration. 
DO  is  lost  through  algal  respiration,  nitrification,  and  breakdown  of 
organic  carbon  (in  particular,  DOC),  and  bottom  sediments.  Five  state 
variables  are  included  in  the  DO  balance:  ammonium,  nitrate,  DOC, 
phytoplankton,  and  DO.  All  constituents  exerting  an  oxygen  demand  must 
be  included  in  the  DO  kinetic  equation. 


The  kinetic  rate  equation  for  in-stream  DO  concentration  is: 


dDO 

dt 


ch 


K{DOs -DOch)  +  roa  (rocaPap  +  rocn(l  -  Pap))ppAp  —  roarocakrpAp 

y 

F7—  {roaVbA  -roakrbAb)-ronFOXnaknitfdTNH+Ach  - rocaFoxmckmcDOCch  -SS0D 

nch 

(53) 


where: 

DOch  =  concentration  of  in-stream  DO  [M/L3] 

ka  =  temperature-dependent  oxygen  reaeration  coefficient  [T1] 
DOs  =  saturation  concentration  of  oxygen  [mg02/L] 

Ssod  =  sediment  oxygen  demand  rate  [M/Ls], 
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Oxygen  reaeration  coefficient 

Oxygen-deficient  (i.e.  below  saturation)  waters  are  replenished  via 
atmospheric  reaeration.  The  reaeration  coefficient  is  computed  as  a 
function  of  the  stream’s  hydraulics  and  (optionally)  wind  velocity: 

K=K,+t^  (54) 


where: 

kah  =  reaeration  rate  based  on  the  stream’s  hydraulic  characteristics 
[T1] 

Kaw  =  reaeration  mass-transfer  coefficient  based  on  wind  velocity 
[L/T], 

Various  methods  have  been  used  to  calculate  atmospheric  reaeration 
coefficients  and  experience  has  shown  that  the  most  effective  method 
depends  upon  the  prevalent  hydraulic  characteristics  of  the  system.  Each 
method  corresponds  to  an  empirical  formula,  which  has  proven  accurate 
for  a  particular  set  of  hydraulic  conditions. 


The  hydraulic  reaeration  equations  pertaining  to  rivers  are  shown  in 
Table  5. 

The  wind  effect  reaeration  equations  are  listed  in  Table  6. 
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Table  5.  Hydraulic  reaeration  equations 


No. 

Equation 

Applicability 

Reference 

1 

°.s 

fc„»=3.93ft,s 

u  -  mean  water  velocity  [m/s] 
h  =  mean  water  depth  [m] 

depths  between 

1-30  ft  and 
velocities 

between  0.5- 1.6 
fps 

O’Connor  and 
Dobbins  (1958) 

2 

k,h  =  5.026  7 

depths  between 

2-11  ft  and 
velocities 

between  1.8-5  fps 

Churchill  et  al. 
(1962) 

3 

Q.e? 

**=5-32^ 

depths  between 

0.4-2. 4  ft  and 

velocities  between 

0.1-1. 8  fps 

Owens  et  al. 

(1964) 

4 

kah=88(u-s)0313h-0353 
for  Q  <  0.556  cms  (<  19.64  cfs) 
kah  =  lA2(u-s)0333h-°-66B;0-243 
for  Q  >  0.556  cms  (>  19.64  cfs) 

for  channel  control 

streams 

Melching  and 

Flores  (1999) 

5 

kah=  2.16(1  +  9Fd0'25)^- 

nch 

U ,  =  JgRhS  is  shear  velocity  [m/s] 

Fd  =  uf \]ghd  is  Froude  number 

Rh  =  hydraulic  radius  [m] 

s  =  channel  slope 

hd  =  the  hydraulic  depth  [m] 

Bt=  top  width  of  the  channel  [m] 

Thackston  and 
Dawson  (2001) 

Table  6.  Wind  effect  reaeration  equations 


No. 

Equation 

Reference 

1 

Kaw  =  0.7281^—0.3171^  +0.0372<10 

uw, io  =  wind  speed  measured  10  meters  above  the  water 
surface  [m/s] 

Banks  and  Herrera  (1977) 

2 

Kaw=  0.09861^ 

Wanninkhof  et  al.  (1991) 

Oxygen  attenuation  coefficient 

For  an  individual  nutrient,  two  equations  are  used  to  represent  the  oxygen 
attenuation  of  the  oxidation  rate: 
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Saturating  exponential 


Half-saturation 


77  (' \  r>  k  oxmnD°ch\ 

b oxmn  =  U  -  e  ) 

(55a) 

Foxna=t  l-ek°™D°ch) 

(55b) 

Foxdn=^-ekoxdnD°ch) 

(55c) 

Foxmp=a-ekoxmpD°ch ) 

(55d) 

F0xmc  =  (1-e  koxmcD°ch) 

(55e) 

Foxmn  =  F>Och  /(Koxmn  +  DO ch) 

(56a) 

Foxna  =  F>Och  / [Koxna  +  DOch  ) 

(56b) 

Foxdn  =  DOch  j (Koxdn  +  DO ch ) 

(56c) 

FOXmp  =  F>Och ! (Koxmp  -)-  DOch ) 

(56d) 

FOXmc  =  F>Och  /  (Koxmc  +  DOch ) 

(56e) 

where 


l^oxmri 

ksxna 

koxdn 

koxmn 

koxmc 

Koxmn 

Koxna 

Koxdn 

Koxmp 

Koxmc 


exponential  coefficient  for  DON  mineralization  [L/mg02] 
exponential  coefficient  for  nitrification  [L/mg02] 
exponential  coefficient  for  denitrification  [L/mg02] 
exponential  coefficient  for  DON  mineralization  [L/mg02] 
exponential  coefficient  for  DOC  mineralization  [L/mg02] 

DO  half-saturation  constant  for  DON  mineralization  [mg02/L] 
DO  half-saturation  constant  for  nitrification  [mg02/L] 

DO  half-saturation  constant  for  denitrification  [mg02/L] 

DO  half-saturation  constant  for  DOP  mineralization  [mg02/L] 
DO  half-saturation  constant  for  DOC  mineralization 
[mg02/L]. 
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Sediment  oxygen  demand 

Sediment  oxygen  demand  (SOD)  is  a  rate  at  which  oxygen  is  consumed  by 
sediments  from  the  overlying  water  column.  It  is  well  known  that  SOD  is 
the  largest  dissolved  oxygen  (DO)  sink  within  natural  waters.  In  addition 
to  chemical  and  biological  processes  within  the  sediment  (DiToro  2001), 
the  rate  of  SOD  also  depends  on  the  water  velocity  above  the 
sediment/water  interface  (Boudreau  and  Joergensen  2001).  Both 
theoretical  (Nakamura  and  Stefan  1994)  and  experimental  (House  2003) 
results  have  provided  relationships  that  describe  the  effect  of  water 
velocity  on  SOD.  In  current  NSM,  SOD  is  modeled  with  a  single 
parameter.  This  simplification  will  be  updated  in  the  future  research. 

Processes  responsible  for  SOD  are  far  too  complex  to  be  represented  by  a 
single  parameter.  In  history,  Fillos  and  Molof  (1972)  modeled  the  SOD 
with  monod  kinetics,  in  which  SOD  is  limited  by  the  DO  in  the  water 
column.  Walker  and  Snodgrass  (1985)  presented  a  slightly  more  complex 
SOD  expression  in  which  SOD  is  divided  into  two  essential  components, 
biological  (BSOD)  and  chemical  (CSOD).  The  former  is  described  by 
monod  kinetics;  the  latter  is  accounted  for  by  first-order  kinetics. 
Higashino  and  Stefan  (2005)  formulated  and  solved  a  numerically 
unsteady  model,  which  relates  SOD  to  sediment  processes  and  bulk  water 
oxygen  by  coupling  diffusive  transport  of  oxygen  in  two  adjacent  layers  at 
the  sediment-water  interface.  In  their  model,  transfer  of  oxygen  to  the 
sediments  is  represented  by  diffusion  through  a  boundary  layer  separating 
the  interface  from  bulk  water,  whereas  penetration  into  the  sediment  layer 
is  accounted  for  by  diffusion  with  oxygen  uptake  represented  by  monod  or 
zero-order  kinetics.  Using  steady-state  compartmental  models,  Di  Toro 
(2001)  developed  steady-state  solutions  for  methane  and  nitrogen,  and 
chemical  SOD,  by  solving  a  coupled  system  of  diffusive-reactive  equations 
expressed  in  terms  of  the  sediment  depth.  Di  Toro  (2001)  provides  a 
comprehensive  analysis  of  primarily  compartmental  dynamic  sediment 
flux  models  for,  among  others,  nutrients  and  sulfide  in  a  sediment  layer 
with  distinct  redox  zones.  Recently,  Hantush  (2007)  derived  unsteady 
solutions  for  pore-water  nitrogen  and  methane  and  their  dissolved  and 
gaseous  fluxes,  and  oxygen  consumption  by  oxidation  reactions  in  sed¬ 
iments  receiving  flux  of  settling  organic  matter  and  interacting  with  the 
bulk  water  through  a  diffusional  boundary  layer.  These  physically  based 
computation  algorithms  will  be  included  in  the  future  version. 
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Appendix  A:  One-Dimensional  Water  Quality 
Transport  Equations 


For  water  quality  modeling  in  streams,  the  advection-dispersion  equation 
with  source/sink  terms  is  one  of  the  governing  mass  balance  equations  to 
be  solved  numerically.  The  primary  purpose  of  this  appendix  is  to  derive 
the  one-dimensional  (lD)  governing  equations  presented  in  the  report. 

Any  water  quality  constituent  in  the  stream  is  transported  by  the  water 
flow  (advection  processes)  and  its  concentration  is  altered  by  the  simul¬ 
taneous  influence  of  turbulent  diffusion  processes.  The  sedimentation  of 
contaminated  suspended  sediments  and  the  bottom  erosion  are  also 
important  pathways  of  the  “water  column-bottom”  exchange  of  constit¬ 
uents.  The  foundation  of  the  channel  water  quality  transport  model  is  the 
lD  mass  balance  equation.  The  mass  balance  for  each  state  variable  tracks 
all  sources,  losses,  and  internal  transformations  of  water  quality  constit¬ 
uents  in  the  stream.  For  each  water  flow  control  volume  (Figure  Ai),  the 
geometry  is  specified  by  a  cell  width  B  length  Ax  and  thickness  h. 


lD  mass  conservation  implies  that: 


d(VCk) 

dt 


-  Advection  +  Dispersion  +  Lateral  Flow  +  Erosion/ Deposition  ^  j_) 
+  External  Sources / Sinks  +  Internal  Sources / Sinks  [Kinetics) 


where: 

Ck  =  concentration  of  constituent  k  in  the  channel  water  [M/Ls] 
V  =  volume  of  the  main  channel  segment  [Ls], 
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Figure  Al.  Water  quality  constituent  flux  through  a  ID  control  volume 

Cross-sectionally  averaged  variables  are  often  used  in  channel  models.  Key 
processes  of  Equation  Al  that  may  be  part  of  any  water  quality  transport 
models  are: 

•  advection-dispersion  transport  of  water  flow 

•  sedimentation  and  resuspension 

•  kinetic  processes  and  external  sources/sinks 

The  net  mass  advected  into  the  control  volume  and  out  of  the  control 
volume  is  given  by: 


„  J  .  d(QCJ 

Advection  = - —Ax 

dx 


(A2) 


where: 

Q  =  volumetric  flow  rate  [L3/T] 
Ax  =  lateral  distance  [L]. 


The  net  dispersive  flux  across  the  control  volume  is: 


Dispersion 


d_ 

dx 


AD,  d°k 


dx 


Ax 


(A3) 


where: 

A  =  channel  cross  section  area  of  the  flow  [L2] 

Dx  =  dispersion  coefficient  of  constituent  j  in  the  x-direction  [L2/T]. 
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The  net  mass  entering  the  channel  via  lateral  inflow  and  leaving  through 
lateral  outflow  is: 


Lateral  Flow  =  ( qlinCkl  -  qloutCk )  Ax  (A4) 


where: 

Cki  =  concentration  of  constituent  k  in  the  lateral  inflow  [M/L'fl 
qiin  =  lateral  inflow  rate  [L2/T] 
qiout  =  lateral  outflow  rate  [L2/T]. 


The  net  upward  flux  of  particulate  substances  near  the  bed  is  the 
difference  between  erosion  and  deposition: 


Erosion/  Deposition  =  A xB(Sr-S 


(  N 


=  A xB 


N 


Y,fpn2VrCk2-Y,fpnVsCk 


(A5) 


where: 

Ck2  =  concentration  of  constituent  k  in  the  upper  layer  of  the  bed 
sediment  [M/L.3] 

Sr  =  rate  per  unit  bed  area  at  which  particulate  substance  is 
entrained  from  the  bed  into  the  water  column  [M/L2/T] 

Ss  =  rate  per  unit  bed  area  at  which  particulate  substance  is 
deposited  from  the  water  column  onto  the  sediment  bed 
[M/L2/T] 

B  =  width  of  the  channel  [L] 

vr  =  erosion  (resuspension)  velocity  of  particulate  substances  [L/T] 
vs  =  effective  settling  (deposition)  velocity  of  particulate  substances 
[L/T], 


Apart  from  the  transport  and  sedimentation  processes,  biochemical 
transformations  significantly  affect  the  concentration  distribution  of  the 
total  substance  in  the  stream.  Major  transformation  processes  that  should 
be  considered  in  a  stream  network  system  have  been  specified  for  each 
water  quality  state  variable  in  NSM.  Substitution  of  all  terms  into 
Equation  At  and  dividing  by  Ax  leads  to  the  following  lD  transport 
equation  for  each  dissolved  state  variable: 
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d(ACk )  _ 

d(QCk)  d 

dt 

dx  dx 

AD 


+  {Qlin^kl  ~  Qlout^k)  +  +^^k,m 


(A6) 


where: 

Rk,m  =  source/sink  term  representing  internal  gains  and  losses  of 
constituent  k  due  to  kinetic  reaction  m  [M/L2/T] 

Sk,m  =  source/sink  term  representing  external  gains  and  losses  of 
constituent  k  [M/L2/T]. 


If  erosion  and  deposition  of  suspended  sediment  occur,  the  following 
equation  for  each  particulate  state  variable  is  appropriate: 


dt 


d(QCk) 

hA[ 

dx 

dx 

AD , 


dc. 


dx 


~b  [QlinP kl  Qlout^k 


+  B 


N  N 

EfPn2VrCk2 -EfpnVsCk  +E  Sk,m+^Rk„ 


(A7) 


Flow  in  the  channel  network  is  modeled  by  the  lD  Saint-Venant  equations. 
The  continuity  equation  (water  balance  law)  accounts  for  the  volume  of 
water  in  a  reach  of  open  channel  including  inflow,  outflow,  and  storage 
within  the  reach. 


The  lD  continuity  equation  is: 


dA  dQ  , 

-^7  =  -  —  +  (Ql,n-Qloul) 

at  ox 


(A8) 


Taking  into  account  the  continuity  equation  and  Equations  A2  and  A3: 


1 

-Jid 

O 

QdCk 

1  d 

dt 

A  dx 

^  A  dx 

AD,  °Ck 


dx 


(A9) 


+j(c»-el)+Es.,+!:s 


k,m 
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+  T  2_^/Jpn2Vr(Jk2  ~  2_^SpnVsCk  +  2_^^k,m  +  ^Rk,1 
“1  1  1 


(A10) 


The  general  mass  balance  equation  for  lD  transport  considering 
advection-dispersion  with  source/sink  terms  is  rewritten  as: 


5C_  QdC  1_ d_ 
dt  A  dx  A  dx 


(All) 


where: 


y^Sk  and  J2Rk  denote  total  external  and  internal  (kinetic)  source/sink 
terms,  respectively. 

In  order  to  provide  a  complete  description  of  the  dominant  pools  and 
fluxes  in  the  water  column,  NSM  takes  into  account  water  quality 
transport  processes  in  the  bed  sediment  (Figure  A2).  The  mass  balance 
approach  is  greatly  simplified  as  the  processes  of  advection,  dispersion, 
and  lateral  inflow  do  not  affect  the  bed  sediment  layer.  For  each  upper 
layer  of  bed  sediment,  water  quality  dynamics  in  the  upper  active  layer  of 
bed  sediments  for  each  dissolved  state  variable  and  each  particulate  state 
variable  can  be  described  as: 


(A12) 


(A13) 


Figure  A2.  Schematic  representation  of  benthic  sediment 
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Appendix  B:  One-Dimensional  Finite 
Difference  Equations  and  Numerical 
Solutions 

The  transport  equations  described  in  Appendix  A  must  be  solved 
simultaneously  for  all  phases  of  the  constituent.  It  is  not  possible  to  solve 
the  coupled  equations  analytically  for  a  general  case.  Therefore  numerical 
techniques  must  be  used.  This  appendix  is  devoted  to  the  formulation  of 
finite  difference  (FD)  equations  and  numerical  solutions.  The  FD 
equations  are  formulated  from  the  original  partial  differential  equations 
(PDEs). 

Finite  difference  approximation  equations 

The  general  mass  balance  equation  of  lD  transport  considering  advection- 
dispersion  with  source/sink  terms  for  a  dissolved  or  suspended  substance 
represented  by  the  concentration  C  is: 

dc_ 

dt 

where: 

C  = 

Ux  = 

^Sk  and  T,Rk  = 

A  common  method  of  solving  PDEs  is  to  approximate  the  spatial 
derivatives  ( of  ox  and d/dy )  using  finite  difference  methods  (FDMs).  The 

FDM  can  be  used  to  solve  these  equations  numerically  in  a  relatively 
simple,  convenient,  and  economical  manner.  To  implement  an  FD  scheme, 
the  channel  reach  is  first  spatially  discretized  into  a  computational  domain 
composed  of  a  number  of  non-overlapping  control  volumes  within  which 
mass  is  conserved.  Figure  Bi  depicts  a  discretized  stream  network  system. 


dc  id 

-u„  —  + 


■'  dx  a  dx 


AD, 


dC 

dx 


+y 


(Bl) 


concentration  (dissolved  or  particulate  phase)  [M/L3] 
depth  averaged  x-direction  (longitudinal)  velocity  [L/T] 

external  and  internal  source/sink  terms,  respectively. 
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Figure  Bl.  Conceptual  stream  network  system 


Referring  to  Figure  B2,  the  first  derivatives  representing  the  advection  and 
dispersion  terms  of  Equation  Bl  for  cell  (j)  can  be  approximated  as 
follows: 


Advection  = 


—u 


dC_ 

dx 


=  ~u« 


Cj+l/2  Cj  1/2 

Ax, 


(B2) 


Dispersion  = 


A  dx 


AD  — 
dx 


1 

x&eJ 

J+l/2 

1  axj 

(B3) 


j- 1/2 


A 


Ax,. 


In  order  to  outline  the  derivation  of  the  FDM,  one  must  start  with  theiD 
advection-dispersion  equation. 


The  FD  equation  for  lD  transport  considering  advection-dispersion 
without  the  source/sink  term  can  be  written  as: 


sin+l  sin  si*  si* 

S’  ~  S'  _  *  ^j+l/2  ~Ui-l/2 

dt  xj  Ax. 


adA 

f  ;/r  1 

dx  , 

j+ 1/2 

(B4) 


I-1/2 


AJ 


AXj 
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where: 


j 

j- 1  andj+i 
j-V 2  andj+V 2 

■5f 


current  location 

positions  removed  and  increased  a  value  x  from  the 
current  location  in  the  solution  scheme,  respectively 
left  and  right  interfaces  of  cell  (j),  respectively 
appropriate  time  level,  which  will  be  determined  later. 


j-1  j 

j+1 

Ik.  v 

j-1/ 2  j+1/2 

1  A*,-i  |  Ax,  ,  Ax,+1 

1  1  1 

1 

1  ^  X 

Figure  B2.  ID  FD  computational  mesh 


Equation  B4  is  the  standard  form  of  the  FD  approximation.  The  interface 
concentration  values  depend  upon  the  appropriate  interpolation  pro¬ 
cedure  for  the  FD  scheme.  How  the  interface  concentration  is  determined 
is  what  distinguishes  one  solution  technique  from  another.  In  the  standard 
finite-difference  method,  the  interface  concentration  is  evaluated  using 
either  the  upstream  or  the  central-in-space  weighting  scheme.  For  the 
upstream  weighting  scheme,  the  interface  concentration  between  two 
neighboring  nodes  is  set  equal  to  the  concentration  at  the  upstream  node 
along  the  same  direction: 


ifuxJ  1/2  >  0 

Cj,  tfuxj_  1/2<0 


(B5) 


Most  FDMs  for  calculating  the  advection  term  are  plagued  by  problems  of 
numerical  oscillations  and/or  artificial  or  numerical  dispersion.  The 
upstream  weighting  scheme  results  in  oscillation-free  solutions.  However, 
the  solution  of  the  advection  term  is  only  accurate  to  the  first  order,  which 
can  lead  to  significant  numerical  dispersion,  since  the  truncation  error 
resulting  from  the  advection  solution  is  of  the  same  order  and  could 
overwhelm  the  second-derivative  physical  dispersion  term  (Zheng  and 
Bennett  1995). 

For  the  central-in-space  weighting  scheme,  the  interface  concentration  is 
set  equal  to  the  weighted  average  of  the  concentrations  on  the  two  sides  of 
the  interface: 
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C, 


i- 1/2 


^  c,  .  +  AXj~1.  C,  =  fl-a.„;,)C,_1+gw_i,2C,(B6) 


Axj._1  +  Ax,.  7  1  Ax,.  .  Ax,.  J 


/0  —  .  J .  C;  +  .  J  C-+1  — (l  a  .+i/2)a+a  .+1/2C,.+1(B7) 


0+1/2 


Ax.+Ax.+1  ^  '  Ax,.+Ax,.+1  J+1 


dc 


<9x 


C,-C/-i 


j_  1/2  0.5(AXj_1  +  AXj) 


(B8) 


3C 


<9x 


C/++-C, 


j+i/2  0.5(Ax,-  +  Ax,-+1) 


(B9) 


where: 


°h/+l/2 


Ax, 


Ax,.+Ax,.+1 


(BIO) 


axj-l/2 


Ax,._± 


Ax  .  ±  +  Ax  . 


(Bll) 


Substitution  of  all  terms  into  Equations  B2  and  B3  leads  to  the  following 
equations: 


Advection 

1 

-u  ■ 


Ax,. 


Ax,.,.  Ax, 

J+1  -C..+ - J- - C 


Ax,. 


Ax,.+AxJ+i  '  Ax,.+AxJ+1 


j+ 1 


-C/-1 


Ax,  1 

-C; 


Ax,-_1  +  Ax,-  J  Ax/_1  +  Ax;.  J 


1  / 


■^4(1-^+i/2)^  +a,i+1/2C,.+1  -(l-a^jC,..,  -a^_1/2C. 


(B12) 


Dispersion 


1  1 


Ax,  Aj 


ADx)i 


Cj+1  Cj  ■( ADJ ,,  C/  Cj  1 


J+1/2  0.5(Ax;. +1  +  Ax;)  v  2  1/2  0.5(Ax;. +  Ax;.) 


j- 1  j  ' 

(B13) 
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Time  weighting  and  central-in-space  weighting  scheme 

Substituting  the  FD  approximations  into  the  PDEs  yields  a  set  of  algebraic 
equations.  Depending  upon  the  manner  in  which  the  FDs  are  computed, 
the  algebraic  equations  may  be  solved  with  either  an  explicit  or  an  implicit 
scheme.  With  an  explicit  scheme,  the  unknown  values  are  found 
recursively.  The  results  of  one  computation  are  necessary  for  the  next. 
With  an  implicit  scheme,  all  the  unknown  values  for  a  given  time  are 
found  simultaneously.  There  are  limitations  for  the  explicit  method  based 
on  the  stability  criterion  associated  with  them,  so  the  size  of  time  steps 
cannot  exceed  a  certain  value.  Due  to  these  limitations,  the  explicit  scheme 
is  not  commonly  used  in  constituent  transport.  As  a  result,  one  would 
naturally  consider  their  implicit  formulations  as  an  alternative  to  avoid 
unstable  solutions.  The  implicit  space-centered  finite  difference  scheme  is 
second-order  accurate  in  both  space  and  time  with  a  small  truncation 
error.  With  the  central -in-space  weighting  scheme  the  solution  of  the 
advection  term  is  accurate  to  the  second  order  and  as  a  result  it  does  not 
lead  to  any  numerical  dispersion  since  the  solution  of  the  dispersion  term 
is  also  accurate  to  the  second  order.  However  if  the  transport  problem  is 
advection-dominated,  the  central-in-space  weighting  scheme  can  lead  to 
excessive  artificial  oscillation,  which  is  typical  of  higher-order  truncation 
errors.  Pinder  and  Gray  (1977)  noted  that  numerical  solutions  of  the 
advection-dispersion  equation  are  characterized  by  two  principal 
phenomena:  numerical  dispersion  and  oscillation  (over-  or  undershoot) 
with  these  phenomena  being  closely  related.  When  a  numerical  scheme  is 
developed  to  minimize  the  numerical  dispersion,  oscillation  is 
encountered  but  when  the  oscillation  is  controlled,  it  is  generally  at  the 
expense  of  increased  numerical  dispersion.  Numerical  dispersion  and 
artificial  oscillation  limit  the  reliability  of  numerical  solutions  of  the 
advection-dispersion  equation  when  FDMs  are  used  (Wang  and  Lacroix 
1997).  Because  of  the  dual  problems  of  numerical  dispersion  and  artificial 
oscillation,  the  standard  FDMs  are  suitable  only  for  solving  transport 
models  not  dominated  by  advection  (i.e.,  when  either  the  physical 
dispersivity  is  large  or  the  grid  spacing  is  made  sufficiently  fine). 

Applying  time  weighting  and  a  central-in-space  weighting  scheme,  the 
spatial  derivatives  are  made  at  both  the  present  and  the  future  time  step. 
These  FD  approximations  are  then  averaged  to  obtain  a  spatial 
approximation  that  corresponds  to  the  midpoint  of  the  time  step. 

Equation  B4  at  any  finite  difference  cell  (/)  can  be  expressed  as  follows: 
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where: 


co  =  time  weighting  coefficient  for  semi-implicit  schemes 
(0  <  co  <  1). 

When  co  =  o  the  solution  scheme  is  explicit  and  when  co  =  l  the  scheme  is 
implicit.  If  =  0.5 ,  Equation  B14  becomes  the  well-known  Crank- 
Nicolson  (C-N)  FD  scheme  approximation. 

Substitution  of  all  interface  terms  into  Equation  B14  leads  to  the  following 
equation: 


c; 


n+ 1 


-c; 


At 


(l-a*+1  /2)q+1 

-(1-a  ■  1/2)C 


_1_  ( y  /nrrt+1 

^  Uxj+l/2^j+l 


n+1 


'  axj-l/2^j 


n+1 


—  (!  —  <*>) 


Ax,. 


- u : 


axj+l  12)^ j  +  axj  1/2C./  1 1 
_ —  axj  1/2  )^i-l  ~  a^7  1/2^- 


+  CO 


a Xj  a;+1 


(ADA 


n+1 


sin+ 1 

S'+l 


-c; 


n+1 


xJj+1/2 


-(ADjr1 


OAfAXj.  +  Axj+1) 

nn+l  _ sin+1 

W  ^J-l 


xJj-1/2 


0.5(Ax-  ±  +  Ax.) 


A-  (1  —  co) 


1  1 


(ADx)%/2 


j-i 

/nrn  /nrn 

^i+i  “S' 


Ax,  A"  , 

J  J  — (AD.)n+1 


0.5(Ax;.  +Axi+1) 

/nrn  _ /nrn 

S'  _S-i 


'xJj-l/2 


OAfAx^  +  Ax,.) 


(B15) 


This,  in  turn,  may  be  simplified  by  collecting  terms: 
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where: 

F,  F,  and  G  are  matrix  coefficients  and  55  is  a  forcing  function. 

For  the  interior  segments  of  the  channel  (from  j  =  2  to  j  =  m-i)  the  matrix 
coefficients  ( E ,  F,  G)  and  the  forcing  function  (55)  are  computed  as 
follows: 
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j+i 


The  application  of  the  time  weighting  and  central-in-space  weighting  finite 
difference  schemes  for  all  inner  cells,  within  the  watershed,  results  in  a  set 
of  tridiagonal  linear  equations.  These  equations  must  be  solved 
simultaneously  to  obtain  the  concentration  C"  1 . 


A  set  of  tridiagonal  equations  representing  M  cells  system  is  shown  below: 
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The  source/sink  terms  must  be  specified  for  each  water  quality  component 
in  order  to  solve  transport  Equation  Bi. 

Boundary  conditions 

Boundary  conditions  must  be  defined  at  upstream  and  downstream 
boundary  cells.  Through  these  cells,  material  is  exchanged  with  the 
environment  outside  the  channel  reach.  An  upstream  boundary  condition 
is  shown  in  Figure  B3,  where  the  interfacial  concentration  at  the  upstream 
boundary  is  set  to  a  user-defined  concentration  Cus. 
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Figure  B3.  Upstream  boundary  condition 


For  the  first  segment  of  the  channel  (j  =  1),  the  matrix  coefficients  (E,  F,  G ) 
and  the  forcing  function  (SS)  are  computed  as  follows: 
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A  downstream  boundary  condition  is  shown  in  Figure  B4  where  the 
dispersive  flux  at  the  downstream  boundary  is  set  to  a  user-defined  flux, 
Fluxds. 
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In  many  applications  the  flux  represented  by  the  downstream  boundary 
condition  is  set  to  zero.  Using  the  definition  of  the  interfacial  first 
derivative  given  by  Equation  B13  and  assuming  Axm+l  =  Axm ,  the 

concentration  in  the  cell  adjoining  the  boundary  Cm+i  is  computed  as: 
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Figure  B4.  Downstream  boundary  condition 

For  the  last  segment  of  the  channel  (j  =  m)  the  matrix  coefficients  (E,  F,  G) 
and  the  forcing  function  (SS)  are  computed  as  follows: 
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Solution  of  Equation  Bi6  requires  that  the  flow  area  and  the  cross- 
sectional  average  velocities  for  each  node  be  obtained  first  from  the 
solution  of  the  flow  equations.  If  the  flow  model  itself  does  not  provide 
flow  velocities  directly,  they  must  be  derived  by  dividing  the  flow  rates  Qj 
by  the  flow  areas  A;  for  each  node,  e.g.,  uxj  =  Qj/Aj.  The  procedure  applied 
at  each  time  step  is  to  first  solve  the  flow  equations  and  then,  given  the 
flow  areas  and  velocities,  solve  the  advection-dispersion  transport 
equation  to  obtain  concentrations. 

Thomas  algorithm 

Tridiagonal  coefficient  matrix  Equation  Bi6,  subject  to  appropriate  initial 
and  boundary  conditions,  can  be  solved  using  an  efficient  formulation  of 
the  Thomas  algorithm  (Chapra  and  Canale  2002).  The  Thomas  algorithm 
is  a  very  efficient  way  of  solving  tridiagonal  linear  systems  and  is  therefore 
the  solver  of  choice  for  most  lD  models.  It  is  also  useful  for  2D  models  that 
have  been  discretized  by  a  finite  difference  scheme  leading  to  a  tridiagonal 
linear  system.  The  matrix  itself  is  not  stored;  only  three  vectors  {E,  F,G) 
are  stored.  These  hold  the  values  on  the  matrix  diagonals. 
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For  a  tridiagonal  matrix  equation:  =  {b} ,  the  Thomas  algorithm 

computational  procedure  consists  of  three  steps: 

l.  Decomposing  [A]  into  a  product  of  lower-diagonal  [l]  and  upper- 
diagonal  [U]  matrices:  [L][f/]  =  [A]. 
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By  considering  the  multiplication  of  [L]  and  [U ] : 

k  =  ei/ti-i 

ti  =  fi-rigi- 1 

Thus  the  values  of  r  and  t  may  be  computed  in  the  order  l,  2,...,  n. 

2.  Forward  substitution  is  through  an  intermediate  vector 
{d}  :  [L]{d}  =  {6} ;  thus  {d}  =  [Lp  {b} : 
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Inspecting  the  multiplication  allows  inversion  of  L  : 
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di  =  bi-ridi_  i 

Thus  values  of  {c/}  may  be  computed  in  the  order  l,  2,...,  n. 

3.  Backward  substitution  to  give  the  final  solution:  [t/]{x}  =  {c/} 
thus]*}  =  [U]  1 {d} : 
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Again,  inspecting  the  multiplication  gives  the  effective  inverse  of  [U 


Xn=dn/tl 


Thus  values  of  {*}  maybe  computed  in  the  order  n,  n-i,...,  1. 

Numerical  implementation 

Many  numerical  approaches  could  be  implemented  to  solve  the  coupled 
FD  equations.  All  sources/sinks  for  water  quality  modeling  are  separated 
into  two  groups:  internal  and  external  terms.  The  division  of  terms  allows 
kinetic  sources/sinks  to  be  updated  at  different  frequencies  than  external 
sources/sinks  -  consistent  with  the  coarser  time  scales  associated  with 
biological  and  chemical  processes  as  opposed  to  hydrologic  transport.  This 
separation  also  allows  computation  times  to  be  reduced.  The  source/sink 
term  represents  a  mass  rate  of  change  [M/T]  of  a  water  quality  constit¬ 
uent.  By  separating  the  kinetic  components  (NSM)  and  transport  com¬ 
ponents  (HEC-RAS,  ADH,  GSSHA),  the  complete  water  quality  mass 
balance  transport  equations  are  solved  through  a  step-by-step  procedure: 
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1.  Evaluate  internal  sources/sinks  due  to  biochemical  transformations  and 
water-sediment  fluxes  ( . )  by  an  explicit  method.  That  is,  all  internal 

source/sink  parameters  in  the  FD  equation  are  evaluated  at  time  t  ( n ). 

2.  Combine  the  effects  of  internal  sources/sinks  ( )  and  external 
loading  ( S . ).  This  step  completes  computation  of  all  sources/sinks  that 
are  contained  in  the  array  SSi  in  the  FD  equation. 

3.  Solve  the  FD  equation  and  compute  concentration  at  time  t+  t  (n+i ) 
using  various  numerical  schemes.  This  step  combines  the  effect  of  the  two 
main  transport  mechanisms  (advection  and  dispersion)  with  biochemical 
transformation  reaction. 
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